home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / LIBIP / PSAF.C < prev    next >
C/C++ Source or Header  |  1999-09-11  |  6KB  |  256 lines
  1. /* 
  2.  * acm_io.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * PSAF
  11.  *
  12.  * function to compute power_spectrum and autocorrelation by FFT:
  13.  *   for angular bend function (-->tan mode) 
  14.  *   for radial function (-->rad mode)
  15.  *   for both (-->all);
  16.  *   -->note: normalization is to N2_VERTEX!!
  17.  *
  18.  * evaluate bending energy on basis of spectral densities
  19.  *
  20.  */
  21. #include <stdio.h>
  22. #include <stdlib.h>
  23. #include <math.h>
  24. #include "ip.h"
  25.  
  26.  
  27. #define log2(a) (log(a)/ 0.6931471805599)
  28. #define exp2(a) (pow((double)2, (double)a))
  29. #define    AP_N_LOC    0
  30.  
  31. #define    FETCH_ALL    1
  32. #undef CHECK_INPUT
  33.  
  34.  
  35. /*
  36.  * psaf()
  37.  *   DESCRIPTION:
  38.  *     evaluate power spectrum and autocorrelation function
  39.  *   ARGUMENTS:
  40.  *     n2: input vector size (long)
  41.  *     inp: input vector (float*)
  42.  *     cf: output vector for correlation function (float*)
  43.  *     ps: output vector for power spectrum (float *)
  44.  *   RETURN VALUE:
  45.  *     none
  46.  */
  47.  
  48. void
  49. psaf (long n2, float *inp, float *cf, float *ps)
  50. {
  51.   int i;
  52.   long n_pow2;
  53.   double temp;
  54.   float *imagp;
  55.   float *realp;
  56.   float inp_max;
  57.  
  58.   /*
  59.    * check for input vector size = power of 2
  60.    */
  61.   n_pow2 = n2;
  62.   temp = log2 ((double) n2);
  63.   if ((temp - (long) temp) != 0.0)
  64.     n_pow2 = (long) exp2 ((double) ((long) temp));
  65.   if (n_pow2 != n2) {
  66.     printf ("psaf: Input vector length (n2 = %d) not power of two!! (required for FFT).\nExiting...\n", n2);
  67.     exit (1);
  68.   }
  69.  
  70.   if ((realp = (float *) calloc (n2, sizeof (float))) == NULL)
  71.       exitmess ("\npasf:...mem allocation for real failed\n", 1);
  72.   if ((imagp = (float *) calloc (n2, sizeof (float))) == NULL)
  73.       exitmess ("\npasf:...mem allocation for imag failed\n", 1);
  74.  
  75.   inp_max = (float) 0.0;
  76.   /*
  77.    * copy the input array to the real array
  78.    */
  79.   for (i = 0; i < n2; i++) {
  80.     *(realp + i) = *(inp + i);
  81.     if (*(inp + i) > inp_max)
  82.       inp_max = *(inp + i);
  83.   }
  84.  
  85.   /*
  86.    * normalize
  87.    */
  88.   for (i = 0; i < n2; i++)
  89.     *(realp + i) /= inp_max;
  90.  
  91.   /*
  92.    * now do the forward 1-D fft
  93.    */
  94.   fft (n2, realp, imagp, -1);
  95.  
  96.   /*
  97.    * calculate the magnitude of the coefficients
  98.    * to get the power spectrum
  99.    */
  100.   for (i = 0; i < n2; i++)
  101.     *(ps + i) = ((*(realp + i) * *(realp + i)) + (*(imagp + i) * *(imagp + i)));
  102.  
  103.   /*
  104.    *       -->note: 
  105.    *          strictly speaking, power spectrum should not 
  106.    *          be divided by n2
  107.    */
  108.   printf ("\n...power spectrum (n = %ld):\n", n2);
  109.   for (i = 0; i < n2; i++) {
  110.     //*(ps + i) = *(ps + i)/(float)n2;
  111.     printf ("%7.4f  ", *(ps + i));
  112.     if ((i + 1) % 8 == 0)
  113.       printf ("\n");
  114.   }
  115.  
  116.   /*
  117.    * if desired, real and imaginary parts of the (complex) FFT may be printed, 
  118.    */
  119.   if (FETCH_ALL == 1) {
  120.     printf ("\n\n...real part of CFFT:\n");
  121.     for (i = 0; i < n2; i++) {
  122.       printf ("%7.4f  ", *(realp + i));
  123.       if ((i + 1) % 8 == 0)
  124.         printf ("\n");
  125.     }
  126.  
  127.     printf ("\n\n...imaginary part of CFFT:\n");
  128.     for (i = 0; i < n2; i++) {
  129.       printf ("%7.4f  ", *(imagp + i));
  130.       if ((i + 1) % 8 == 0)
  131.         printf ("\n");
  132.     }
  133.   }
  134.   printf ("\n...prepare to compute back transform:");
  135.   /*
  136.    * copy the power spectrum array to the real array
  137.    * and zero the imag array
  138.    */
  139.   for (i = 0; i < n2; i++) {
  140.     *(realp + i) = *(ps + i);
  141.     *(imagp + i) = (float) 0.0;
  142.   }
  143.  
  144.   /*
  145.    * now do the reverse 1-D fft
  146.    */
  147.   fft (n2, realp, imagp, -1);
  148.  
  149.   /*
  150.    * copy the real part to the autocorrelation array
  151.    */
  152.   for (i = 0; i < n2; i++)
  153.     *(cf + i) = *(realp + i);
  154.  
  155.   printf ("\n...output of autocorrelation (n = %ld):\n", n2);
  156.   for (i = 0; i < n2; i++) {
  157.     /*
  158.      * Ordinarily one would multiply by n2, but
  159.      * given our definition of fft(), we don't
  160.      */
  161.     //*(cf+i) *= (float)n2;
  162.     printf ("%7.4f  ", *(cf + i));
  163.     if ((i + 1) % 8 == 0)
  164.       printf ("\n");
  165.   }
  166.   free (realp);
  167.   free (imagp);
  168.   return;
  169. }
  170.  
  171.  
  172. /*
  173.  * power spectrum and correlation function
  174.  *
  175.  * 1D discrete autocorrelation function, using FFT
  176.  */
  177.  
  178. /*
  179.  * spec_and_corr()
  180.  *   DESCRIPTION:
  181.  *     1D discrete autocorrelation function, using FFT
  182.  *     wrapper for psaf
  183.  *   ARGUMENTS:
  184.  *     n2: input vector size (long)
  185.  *     inp: input vector (float*)
  186.  *     cf: output vector for correlation function (float*)
  187.  *     ps: output vector for power spectrum (float *)
  188.  *   RETURN VALUE:
  189.  *     none
  190.  */
  191.  
  192. void
  193. spec_and_corr (long n2, float *inp, float *cf, float *ps)
  194. {
  195.  
  196. #ifdef CHECK_INPUT
  197.   int i;
  198.  
  199.   printf ("\n...AP input data:\n");
  200.   printf ("...n2 = %ld\n", n2);
  201.   for (i = 0; i < n2; i++) {
  202.     printf ("%7.4f  ", *(inp + i));
  203.     if ((i + 1) % 8 == 0)
  204.       printf ("\n");
  205.   }
  206. #endif
  207.  
  208.  
  209.  
  210. /*
  211.  * evaluate power spectrum and autocorrelation function by 1D-FFT
  212.  */
  213.   psaf (n2, inp, cf, ps);
  214.  
  215. }
  216.  
  217.  
  218. /*
  219.  * spec_bend_en()
  220.  *   DESCRIPTION:
  221.  *     evaluate bending energy
  222.  *     ref: Young et al., Information & Control 25, 357-370 (1974)
  223.  *     -- note: 
  224.  *       dependent on parametrization of contour and normalization of
  225.  *       shape function;
  226.  *       strictly defined for radial function only;
  227.  *   ARGUMENTS:
  228.  *     nv: input vector size (long)
  229.  *     ps: input vector for power spectrum (float *)
  230.  *     c_len: curvature length (double)
  231.  *   RETURN VALUE:
  232.  *     bending energy (double)
  233.  */
  234. double
  235. spec_bend_en (nv, ps, c_len)
  236.      long nv;
  237.      float *ps;
  238.      double c_len;
  239. {
  240.   int i, i4;
  241.   double om0, ps2;
  242.   double e_bend;
  243.  
  244.   e_bend = 0.0;
  245.   for (i = 0; i < nv; i++) {
  246.     i4 = i * i * i * i;
  247.     ps2 = (*(ps + i)) * (*(ps + i));
  248.     e_bend += (i4 * ps2);
  249.   }
  250.   om0 = 2.0 * PI / c_len;
  251.   e_bend *= om0 * om0 * om0 * om0;
  252.   printf ("\n...bending energy (from spectrum) = %lf\n", e_bend);
  253.  
  254.   return (e_bend);
  255. }
  256.